###############################################################################################################
# Replication for individual-level analysis using data from ALLBUS survey for 
# Charnysh and Schaub, Migration and Social Change: Evidence from post-WWII Displacement in Germany, JOP
# Access to the geo-referenced version of the ALLBUS data, including the municipality identifiers (AGS) 
# needed to relate the individual-level results to changes in religious diversity (and controls), is available 
# through the GESIS Secure Data Center (SDC) in Cologne. 
# We provide the output reported in the manuscript in the log file log_ALLBUS.txt
################################################################################################################

rm(list = ls())  #JK
library(readstata13)
library(readxl)
library(lme4)
library(stargazer)
library(tidyverse)
library(lfe)
library(sandwich)
library(corrplot)
library(psych)
library(ggplot2)
library(arm) 
library(tidyverse)
library(knitr)
library(dotwhisker)
library(sdcR) #JK

######################################### Exchange for sdcR-package (JK) ######################################### 
setwd("C:/Users/ra000022/ownCloud")
my_log<-file("Output/apr2025_log.txt") #File path may need to be changed
sink(my_log, append=TRUE, type="output")
sink(my_log, append=TRUE, type="message")
cat(readChar(rstudioapi::getSourceEditorContext()$path, file.info(rstudioapi::getSourceEditorContext()$path)$size))
 
######################################### JK ######################################### 


###################################
## READING IN AND MERGING THE DATA:
###################################

dat<-read_excel("Input/User/municipality-dataset-short.xlsx") #data with historical characteristics of Bavarian municipalities in 1939 and 1950 based on census
allbus<-read.dta13("Input/GESIS/Charnysh_5276.dta") #Reading in ALLBUS survey
bavaria<-subset(allbus, land=="BAVARIA" & year>1992) #Subsetting to Bavaria 

bavaria$agsys<-as.numeric(bavaria$agsys) #unique AGS used for the merge
bav9418<-merge(bavaria, dat, by.x="agsys", by.y="ags") #merging municipality dataset to allbus by ags

#computing historical variables (same as in municipality dataset)

bav9418$sprotestants1939<-bav9418$protestants1939/bav9418$pop1939orig
bav9418$scatholics1939<-bav9418$catholics1939/bav9418$pop1939orig
bav9418$sotherrel1939<-(bav9418$pop1939orig - (bav9418$catholics1939 + bav9418$protestants1939))/bav9418$pop1939orig
bav9418$relfrac1939<-1-(bav9418$sprotestants1939^2+bav9418$scatholics1939^2+bav9418$sotherrel1939^2)

bav9418$sprotestants1950<-bav9418$protestants1950/bav9418$pop1950
bav9418$scatholics1950<-bav9418$catholics1950/bav9418$pop1950
bav9418$sotherrel1950<-(bav9418$pop1950 - (bav9418$catholics1950 + bav9418$protestants1950))/bav9418$pop1950
bav9418$relfrac1950<-1-(bav9418$sprotestants1950^2+bav9418$scatholics1950^2+bav9418$sotherrel1950^2)

bav9418$DeltaFrac<-bav9418$relfrac1950-bav9418$relfrac1939
bav9418$lnpop1939<-log(bav9418$pop1939) 
bav9418$srefugees1950<-bav9418$prefugees1950/100 


####################################################
#### COMPUTING KEY VARIABLES FROM ALLBUS QUESTIONS: 
####################################################

bav9418$female<-ifelse(bav9418$sex=="MALE", 0, 1)

#create year-respid variable for merging in results from principal component analysis
bav9418$RowNM <- interaction(bav9418$year,bav9418$respid, sep = "-")
 
##Traditional gender roles:
bav9418$wom1<-ifelse(bav9418$fr01=="COMPLETELY DISAGREE", 1, ifelse(bav9418$fr01=="TEND TO DISAGREE", 2, ifelse(bav9418$fr01=="TEND TO AGREE", 3, ifelse(bav9418$fr01=="COMPLETELY AGREE", 4, NA))))
bav9418$wom2<-ifelse(bav9418$fr02=="COMPLETELY DISAGREE", 4, ifelse(bav9418$fr02=="TEND TO DISAGREE", 3, ifelse(bav9418$fr02=="TEND TO AGREE", 2, ifelse(bav9418$fr02=="COMPLETELY AGREE", 1, NA))))
bav9418$wom3<-ifelse(bav9418$fr03a=="COMPLETELY DISAGREE", 4, ifelse(bav9418$fr03a=="TEND TO DISAGREE", 3, ifelse(bav9418$fr03a=="TEND TO AGREE", 2, ifelse(bav9418$fr03a=="COMPLETELY AGREE", 1, NA))))
bav9418$wom4<-ifelse(bav9418$fr04a=="COMPLETELY DISAGREE", 4, ifelse(bav9418$fr04a=="TEND TO DISAGREE", 3, ifelse(bav9418$fr04a=="TEND TO AGREE", 2, ifelse(bav9418$fr04a=="COMPLETELY AGREE", 1, NA))))
bav9418$wom5<-ifelse(bav9418$fr05a=="COMPLETELY DISAGREE", 1, ifelse(bav9418$fr05a=="TEND TO DISAGREE", 2, ifelse(bav9418$fr05a=="TEND TO AGREE", 3, ifelse(bav9418$fr05a=="COMPLETELY AGREE", 4, NA))))
bav9418$wom6<-ifelse(bav9418$fr06=="COMPLETELY DISAGREE", 4, ifelse(bav9418$fr06=="TEND TO DISAGREE", 3, ifelse(bav9418$fr06=="TEND TO AGREE", 2, ifelse(bav9418$fr06=="COMPLETELY AGREE", 1, NA))))
 
bav9418$WomenAvg6<-rowMeans(bav9418[, c("wom1", "wom2", "wom3", "wom4", "wom5", "wom6")], na.rm=FALSE)
 
#creating index using principal components
women_vars <- bav9418[c("wom1", "wom2", "wom3", "wom4", "wom5", "wom6")]
rownames(women_vars) <- interaction(bav9418$year,bav9418$respid, sep = "-")
colnames(women_vars) <- c("wom1", "wom2", "wom3", "wom4", "wom5", "wom6")
corrplot(cor(women_vars, use = "complete.obs"))
women_pca <- princomp(na.omit(women_vars), cor = TRUE, scores = TRUE)
summary(women_pca)

women_pca_nrot <- principal(na.omit(women_vars), nfactors = 1, rotate = "none")
women_pca_scores <- women_pca_nrot$scores
loadings(women_pca_nrot)

bav9418 <- merge(bav9418, women_pca_scores, by.x="RowNM", by.y = "row.names", all.x = TRUE)
names(bav9418)[which(names(bav9418)=="PC1")]<-"WomenPCA1"
 
##Prohibiting by law:
bav9418$abortion<-ifelse(bav9418$ca16=="DO NOT PROHIBIT", 1,  ifelse(bav9418$ca16=="PROHIBIT", 0, NA))
bav9418$suicide<-ifelse(bav9418$ca17=="DO NOT PROHIBIT", 1,  ifelse(bav9418$ca17=="PROHIBIT", 0, NA))
bav9418$homosex<-ifelse(bav9418$ca20=="DO NOT PROHIBIT", 1,  ifelse(bav9418$ca20=="PROHIBIT", 0, NA))

#Religion and church:
bav9418$churchAttendance<-ifelse(bav9418$rp01=="MORE THAN 1X A WEEK", 6, ifelse(bav9418$rp01=="ONCE A WEEK", 5, ifelse(bav9418$rp01=="1-3 TIMES A MONTH", 4,ifelse(bav9418$rp01=="SEVERAL TIMES A YEAR", 3, ifelse(bav9418$rp01=="LESS OFTEN", 2, ifelse(bav9418$rp01=="NEVER", 1, NA)))) ))

#expellee indicators for ALLBUS respondents:
bav9418$expellee<-ifelse(bav9418$dm01 %in% c("FORMER GERM.TERR.", "SOVIET UNION", "POLAND", "CZECHOSLOVAKIA", "YUGOSLAVIA", "ROMANIA", "HUNGARY") & bav9418$yborn<1950, 1, ifelse(is.na(bav9418$dm01) & bav9418$yborn<1950, NA, 0))
bav9418$expelleeFather<-ifelse(bav9418$fdm01 %in% c("FORMER GERM.TERR.", "SOVIET UNION", "POLAND", "CZECHOSLOVAKIA", "YUGOSLAVIA", "ROMANIA", "HUNGARY") & bav9418$yborn>1945, 1, ifelse(is.na(bav9418$fdm01), NA, 0))
bav9418$expelleeMother<-ifelse(bav9418$mdm01 %in% c("FORMER GERM.TERR.", "SOVIET UNION", "POLAND", "CZECHOSLOVAKIA", "YUGOSLAVIA", "ROMANIA", "HUNGARY") & bav9418$yborn>1945, 1, ifelse(is.na(bav9418$mdm01), NA, 0))
bav9418$expellee_Origin<-ifelse(bav9418$expellee==1 | bav9418$expelleeFather==1 | bav9418$expelleeMother==1, 1, ifelse(bav9418$expellee==0 & bav9418$expelleeFather==0 & bav9418$expelleeMother==0, 0, NA))

##subset for analysis without expellees

bav9418nonexp<-subset(bav9418, expellee==0) #all non-expellee respondents in Bavaria




##################################################################################################################
##### Table A2: Descriptive statistics for outcome variables based on ALLBUS survey responses in Bavaria ########
##################################################################################################################

vars<-c("WomenPCA1", "wom1",  "wom2",  "wom3",  "wom4","wom5",  "wom6", "abortion", "suicide", "homosex", "churchAttendance")

stargazer(bav9418[,which(names(bav9418) %in% vars) ], title="Table A2: Descriptive statistics for Bavarian municipalities in ALLBUS", digits=2, summary.stat=c("n", "mean","sd",  "min", "max"),  label = "tab:descriptive", type="text",
          covariate.labels =c("Working mother (W1)", "Wife's career (W2)", "Small child's mother (W3)", "Woman works (W4)", "Mother's job (W5)", "Married woman's work (W6)", "Principal Component 1 (gender norms)", "Abortion", "Suicide", "Homosexuality", "Frequency of church attendance"))


##################################################################################################################
##### Table 3: Religious diversity and norms on gender, sexuality, doctor-assisted suicide, and abortion #########
##################################################################################################################

gender_mod1r <- felm(WomenPCA1 ~ sex + age +relfrac1950 + relfrac1939 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
gender_mod2r <- felm(WomenPCA1 ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939+srefugees1950 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
homos_felm1r<- felm(homosex ~ sex + age + relfrac1950 +relfrac1939 | 0 |0|agsys, data = bav9418, weights=bav9418$wghtpt)
homos_felm2r<- felm(homosex ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939+srefugees1950 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
suic_felm1r<- felm(suicide ~ sex + age + relfrac1950 +relfrac1939 | 0 |0|agsys, data = bav9418, weights=bav9418$wghtpt)
suic_felm2r<- felm(suicide ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939+srefugees1950 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
abortion_felm1r<- felm(abortion ~ sex + age + relfrac1950 +relfrac1939 | 0 |0|agsys, data = bav9418, weights=bav9418$wghtpt)
abortion_felm2r<- felm(abortion ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939+srefugees1950 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 

stargazer(gender_mod1r, gender_mod2r, homos_felm1r, homos_felm2r, suic_felm1r, suic_felm2r, abortion_felm1r, abortion_felm2r, no.space=TRUE, type="text", 
          title="Table 3: Religious diversity and norms on gender, sexuality, doctor-assisted suicide, and abortion.", digits=2, 
          column.labels =c("minimal","full mun covars", "minimal","full mun covars", "minimal","full mun covars"))


############################################################
### Table A21: Religious diversity and church attendance. ##
############################################################

chur_felm1 <- felm(churchAttendance ~ sex + age +relfrac1950 + relfrac1939 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
chur_felm2 <- felm(churchAttendance ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939| year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
chur_felm2r <- felm(churchAttendance ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939+srefugees1950 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
chur_felm3r <- felm(churchAttendance ~ sex + age +expellee_Origin+relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939+srefugees1950  | year|0|agsys, data = bav9418, weights=bav9418$wghtpt)
chur_felm5r <- felm(churchAttendance ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939+srefugees1950  | year |0|agsys, data = bav9418nonexp, weights=bav9418nonexp$wghtpt) 

stargazer(chur_felm1, chur_felm2, chur_felm2r, chur_felm3r, chur_felm5r, no.space=TRUE, type="text", title="Table: ALLBUS survey. Religious diversity and church attendance.", digits=2,  
          column.labels =c("minimal","full mun covars", "adding expellee control", "adding expellee covars", "no expellees subset" ))


################################################
##### Table A22: Individual gender questions ###
################################################

#Panel A: Main specification

wom_felm1r<- felm(wom1 ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939+srefugees1950  | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
wom_felm2r <- felm(wom2 ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939+srefugees1950  | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
wom_felm3r <- felm(wom3 ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939+srefugees1950  | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
wom_felm4r <- felm(wom4 ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939+srefugees1950  | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
wom_felm5r <- felm(wom5 ~  sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939+srefugees1950  | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
wom_felm6r <- felm(wom6 ~  sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939+srefugees1950  | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 

stargazer(wom_felm1r, wom_felm3r, wom_felm5r, wom_felm2r, wom_felm4r, wom_felm6r, no.space=TRUE, type="text",  title="Table, Panel 1: ALLBUS survey. Religious diversity and gender norms, controlling for expellee share.", digits=2)

#Panel B: Alternative specification without Expellee share.

wom_felm1<- felm(wom1 ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
wom_felm2 <- felm(wom2 ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
wom_felm3 <- felm(wom3 ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
wom_felm4 <- felm(wom4 ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
wom_felm5 <- felm(wom5 ~  sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
wom_felm6 <- felm(wom6 ~  sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 

stargazer(wom_felm1, wom_felm3, wom_felm5, wom_felm2, wom_felm4, wom_felm6, no.space=TRUE, type="text",  title="Table, Panel 2: ALLBUS survey. Religious diversity and gender norms, not controlling for expellee share.", digits=2)


##############################################
## Table A23: Religious diversity and norms ##
##############################################

#basic model and model with municipal covariates, alternative specification without expellee share.
gender_mod1 <- felm(WomenPCA1 ~ sex + age +relfrac1950 + relfrac1939 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
gender_mod2 <- felm(WomenPCA1 ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 

homos_felm1<- felm(homosex ~ sex + age + relfrac1950 +relfrac1939 | 0 |0|agsys, data = bav9418, weights=bav9418$wghtpt)
homos_felm2<- felm(homosex ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
suic_felm1<- felm(suicide ~ sex + age + relfrac1950 +relfrac1939 | 0 |0|agsys, data = bav9418, weights=bav9418$wghtpt)
suic_felm2<- felm(suicide ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 
abortion_felm1<- felm(abortion ~ sex + age + relfrac1950 +relfrac1939 | 0 |0|agsys, data = bav9418, weights=bav9418$wghtpt)
abortion_felm2<- felm(abortion ~ sex + age +relfrac1950 + relfrac1939+sagri1939+sfemale1939+popdens1939+log(distBorderEast)+dest1945+ diststation1950+lnpop1939 | year |0|agsys, data = bav9418, weights=bav9418$wghtpt) 

stargazer(gender_mod1, gender_mod2, homos_felm1, homos_felm2, suic_felm1, suic_felm2, abortion_felm1, abortion_felm2, no.space=TRUE, type="text",  title="Religious diversity and norms on gender, sexuality, doctor-assisted suicide, and
abortion. Alternative specification without expellee share.", digits=2, column.labels =c("minimal","full mun covars", "minimal","full mun covars", "minimal","full mun covars"))


